SRCD 2019 Symposium: Religious & metaphysical concepts (Srinivasan)

# load required libraries
library(tidyverse)

── Attaching packages ─────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.0.0     ✔ purrr   0.2.5
✔ tibble  1.4.2     ✔ dplyr   0.7.6
✔ tidyr   0.8.1     ✔ stringr 1.3.1
✔ readr   1.1.1     ✔ forcats 0.3.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(langcog) # source: https://github.com/langcog/langcog-package

Attaching package: ‘langcog’

The following object is masked from ‘package:base’:

    scale

library(psych)

Attaching package: ‘psych’

The following objects are masked from ‘package:ggplot2’:

    %+%, alpha

library(lme4)

Loading required package: Matrix

Attaching package: ‘Matrix’

The following object is masked from ‘package:tidyr’:

    expand

library(cowplot)

Attaching package: ‘cowplot’

The following object is masked from ‘package:ggplot2’:

    ggsave

# set theme for ggplots
theme_set(theme_bw())

Data preparation

d_all <- # d_us_ad_pilot %>% rownames_to_column("subid") %>%
  d_gh_ad %>% rownames_to_column("subid") %>%
  full_join(d_gh_ch %>% rownames_to_column("subid")) %>%
  full_join(d_th_ad %>% rownames_to_column("subid")) %>%
  full_join(d_th_ch %>% rownames_to_column("subid")) %>%
  full_join(d_vt_ad %>% rownames_to_column("subid")) %>%
  full_join(d_vt_ch %>% rownames_to_column("subid")) %>%
  column_to_rownames("subid")

Joining, by = c("subid", "choose what to do", "feel guilty", "feel happy", "feel love", "feel pain", "feel proud", "feel sad", "feel scared", "feel shy", "feel sick, like when you feel like you might vomit", "feel tired", "figure out how to do things", "get angry", "get hungry", "get hurt feelings", "hear things", "pray", "remember things", "sense temperatures", "sense when things are far away", "smell things", "think about things")
Joining, by = c("subid", "choose what to do", "feel guilty", "feel happy", "feel love", "feel pain", "feel proud", "feel sad", "feel scared", "feel shy", "feel sick, like when you feel like you might vomit", "feel tired", "figure out how to do things", "get angry", "get hungry", "get hurt feelings", "hear things", "pray", "remember things", "sense temperatures", "sense when things are far away", "smell things", "think about things")
Joining, by = c("subid", "choose what to do", "feel guilty", "feel happy", "feel love", "feel pain", "feel proud", "feel sad", "feel scared", "feel shy", "feel sick, like when you feel like you might vomit", "feel tired", "figure out how to do things", "get angry", "get hungry", "get hurt feelings", "hear things", "pray", "remember things", "sense temperatures", "sense when things are far away", "smell things", "think about things")
Joining, by = c("subid", "choose what to do", "feel guilty", "feel happy", "feel love", "feel pain", "feel proud", "feel sad", "feel scared", "feel shy", "feel sick, like when you feel like you might vomit", "feel tired", "figure out how to do things", "get angry", "get hungry", "get hurt feelings", "hear things", "pray", "remember things", "sense temperatures", "sense when things are far away", "smell things", "think about things")
Joining, by = c("subid", "choose what to do", "feel guilty", "feel happy", "feel love", "feel pain", "feel proud", "feel sad", "feel scared", "feel shy", "feel sick, like when you feel like you might vomit", "feel tired", "figure out how to do things", "get angry", "get hungry", "get hurt feelings", "hear things", "pray", "remember things", "sense temperatures", "sense when things are far away", "smell things", "think about things")

Shared conceptual structure

Pooling all participants from all sites together into a common factor structure.

Parallel analysis

How many factors to retain?

reten_all_PA <- fa.parallel(d_all, plot = F); reten_all_PA

Parallel analysis suggests that the number of factors =  4  and the number of components =  2 
Call: fa.parallel(x = d_all, plot = F)
Parallel analysis suggests that the number of factors =  4  and the number of components =  2 

 Eigen Values of 

reten_all_par <- reten_all_PA$nfact

What are these factors?

efa_all_par <- fa(d_all, nfactors = reten_all_par, rotate = "oblimin",
                  scores = "tenBerge", impute = "median")

Loading required namespace: GPArotation

heatmap_fun(efa_all_par, factor_names = c("COGNITIVE", "BODILY", "SOCIAL-EMOTIONAL", "[other]")) + 
  labs(title = "Factor loadings (data pooled across all samples)")

the condition has length > 1 and only the first element will be usedJoining, by = "capacity"
Joining, by = "factor"

Which capacities are attributed to which targets?

Factor scores

scoresplot_fun(efa_all_par, target = "all",
               factor_names = c("COGNITIVE", "BODILY", 
                                "SOCIAL-EMOTIONAL", "[other]")) + 
  labs(title = "Factor scores (by sample, factor, and target entity)",
       subtitle = "All targets")

the condition has length > 1 and only the first element will be usedIgnoring unknown aesthetics: y

scoresplot_fun(efa_all_par, target = c("ghosts", "God", "children"),
               factor_names = c("COGNITIVE", "BODILY",
                                "SOC.-EMO.", "[other]")) + 
  labs(title = "Factor scores (by sample, factor, and target entity)",
       subtitle = "Human and 'supernatural' targets only")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedIgnoring unknown aesthetics: y

itemsplot_fun(efa_all_par, target = c("ghosts", "God", "children")) +
  labs(title = "Parallel Analysis")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedJoining, by = "capacity"

|=============================                          | 53% ~2 s remaining     
|==============================                         | 55% ~2 s remaining     
|===============================                        | 57% ~2 s remaining     
|================================                       | 59% ~2 s remaining     
|=================================                      | 60% ~1 s remaining     
|==================================                     | 62% ~1 s remaining     
|==================================                     | 63% ~1 s remaining     
|====================================                   | 66% ~1 s remaining     
|=====================================                  | 67% ~1 s remaining     
|=====================================                  | 69% ~1 s remaining     
|=======================================                | 71% ~1 s remaining     
|=======================================                | 72% ~1 s remaining     
|========================================               | 74% ~1 s remaining     
|=========================================              | 76% ~1 s remaining     
|==========================================             | 77% ~1 s remaining     
|===========================================            | 79% ~1 s remaining     
|===========================================            | 80% ~1 s remaining     
|============================================           | 81% ~1 s remaining     
|=============================================          | 83% ~1 s remaining     
|==============================================         | 85% ~1 s remaining     
|===============================================        | 87% ~0 s remaining     
|================================================       | 89% ~0 s remaining     
|=================================================      | 91% ~0 s remaining     
|===================================================    | 93% ~0 s remaining     
|====================================================   | 95% ~0 s remaining     
|=====================================================  | 97% ~0 s remaining     
|====================================================== | 99% ~0 s remaining     

Joining, by = c("factor", "capacity", "order")

Summary scores

factors_par <- efa_all_par$loadings[] %>%
  data.frame() %>%
  rownames_to_column("capacity") %>%
  gather(factor, loading, -capacity) %>%
  group_by(capacity) %>%
  top_n(1, loading) %>%
  ungroup() %>%
  arrange(factor, desc(loading)) %>%
  mutate(order = 1:nrow(.))
factors_par_howmany <- factors_par %>%
  count(factor)
factors_par_culled <- factors_par %>%
  group_by(factor) %>%
  top_n(min(factors_par_howmany$n), loading) %>%
  ungroup() %>%
  arrange(factor, desc(loading)) %>%
  mutate(order = 1:nrow(.))
scores_par_prelim <- d_all %>%
  rownames_to_column("site_age_subid_target_entity") %>%
  gather(capacity, response, -site_age_subid_target_entity) %>%
  mutate(site = gsub("_.*$", "", site_age_subid_target_entity),
         age = gsub("_.*$", "", gsub("^.._", "", site_age_subid_target_entity)),
         subid = gsub("_target.*$", "", site_age_subid_target_entity),
         target = gsub("^.*target_", "", site_age_subid_target_entity)) %>%
  mutate(target = case_when(target == "pigs" ~ NA_character_,
                            target == "NA" ~ NA_character_,
                            target == "crickets" ~ "beetles",
                            TRUE ~ target)) %>%
  mutate(site = factor(site, levels = c("us", "gh", "th", "ch", "vt"),
                       labels = c("US", "Ghana", "Thailand", "China", "Vanuatu")),
         age = factor(age, levels = c("ad", "ch"),
                      labels = c("adults", "children")),
         target = factor(target,
                         levels = c("rocks", "flowers", "cell phones", 
                                    "beetles", "chickens", "mice", 
                                    "dogs", "children", 
                                    "ghosts", "God"))) %>%
  select(-site_age_subid_target_entity) %>%
  full_join(factors_par_culled) %>%
  # full_join(factors_par) %>%
  mutate(factor = factor(factor, levels = c("MR2", "MR1", "MR3", "MR4"),
                         labels = c("BODILY", "COGNITIVE", 
                                    "SOCIAL-EMOTIONAL", "[other]"))) %>%
  filter(!is.na(target), !is.na(factor)) %>%
  distinct()

Joining, by = "capacity"

scores_par <- scores_par_prelim %>%
  group_by(site, age, subid, target, factor) %>%
  summarise(score = mean(response, na.rm = T)) %>%
  ungroup()
scores_par_mb <- scores_par %>% 
  group_by(site, age, target, factor) %>%
  multi_boot_standard(col = "score", na.rm = T) %>%
  ungroup()

# calculate reliability
scores_par_prelim_wide <- scores_par_prelim %>%
  full_join(factors_par_culled) %>%
  filter(!is.na(factor), !is.na(subid)) %>% 
  arrange(factor, desc(loading)) %>%
  select(subid, capacity, response) %>%
  spread(capacity, response) %>%
  select(subid, factors_par_culled$capacity) %>%
  column_to_rownames("subid")

Joining, by = c("capacity", "factor", "loading", "order")
Column `factor` joining factor and character vector, coercing into character vector

# BODILY
alpha(scores_par_prelim_wide[5:8])$total

# "COGNITIVE"
alpha(scores_par_prelim_wide[1:4])$total

# "SOCIAL-EMOTIONAL"
alpha(scores_par_prelim_wide[9:12])$total

ggplot(scores_par, aes(x = target, y = score, color = factor)) +
  facet_grid(cols = vars(site, age), rows = vars(factor)) +
  geom_jitter(height = 0.02, width = 0.2, alpha = 0.2, show.legend = F) +
  geom_pointrange(data = scores_par_mb,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  fatten = 2, color = "black") +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

ggplot(scores_par %>% filter(target %in% c("children", "ghosts", "God")), 
       aes(x = target, y = score, color = factor)) +
  facet_grid(cols = vars(site, age), rows = vars(factor)) +
  geom_jitter(height = 0.02, width = 0.2, alpha = 0.25, show.legend = F) +
  geom_pointrange(data = scores_par_mb %>%
                    filter(target %in% c("children", "ghosts", "God")),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  fatten = 2, color = "black") +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

Minimizing BIC

How many factors to retain?

reten_all_BIC <- VSS(d_all, plot = F); reten_all_BIC
reten_all_bic <- data.frame(reten_all_BIC$vss.stats %>% rownames_to_column("nfact") %>% top_n(-1, BIC))$nfact %>% as.numeric()

What are these factors?

efa_all_bic <- fa(d_all, nfactors = reten_all_bic, rotate = "oblimin",
                  scores = "tenBerge", impute = "median")

heatmap_fun(efa_all_bic, 
            factor_names = c("COGNITIVE", "BODILY", "SOCIAL-EMOTIONAL", "other")) + 
  # labs(title = "Minimizing BIC")
  labs(x = "Figure 1: Shared conceptual structure") +
  theme(axis.title.x = element_text(hjust = 0)) +
  guides(fill = guide_colorbar("factor loading", barheight = 15))

Which capacities are attributed to which targets?

scoresplot_fun(efa_all_bic, target = "all", highlight = "supernatural",
            factor_names = c("COGNITIVE", "BODILY", "SOCIAL-EMOTIONAL", "other")) + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 1,
                                   color = c(rep("black", 8),
                                             rep("#984ea3", 2)),
                                   face = c(rep("plain", 8),
                                            rep("bold", 2)))) +
  scale_fill_manual(values = c("#e41a1c", "#4daf4a", "#377eb8", "gray")) +
  # labs(title = "Minimizing BIC")
  labs(title = "Figure 2: Factor scores")

scoresplot_fun(efa_all_bic, target = c("ghosts", "God", "children")) + 
  # scale_fill_manual(values = c("#e41a1c", "#4daf4a", "#377eb8")) +
  labs(title = "Minimizing BIC")

itemsplot_fun(efa_all_bic, target = c("ghosts", "God", "children")) + 
  labs(title = "Minimizing BIC")

Preset criteria

How many factors to retain?

reten_all_k <- reten_fun(d_all, rot_type = "oblimin"); reten_all_k

What are these factors?

efa_all_k <- fa(d_all, nfactors = reten_all_k, rotate = "oblimin",
                  scores = "tenBerge", impute = "median")

heatmap_fun(efa_all_k) + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

Which capacities are attributed to which targets?

scoresplot_fun(efa_all_k, target = "all") + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

scoresplot_fun(efa_all_k, target = c("ghosts", "God", "children")) + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

itemsplot_fun(efa_all_k, target = c("ghosts", "God", "children")) + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

3 factors

What are these factors?

efa_all_3 <- fa(d_all, nfactors = 3, rotate = "oblimin",
                scores = "tenBerge", impute = "median")

heatmap_fun(efa_all_3, 
            factor_names = c("COGNITIVE", "BODILY", "SOCIAL-EMOTIONAL")) + 
  labs(x = "Shared concpetual structure") +
  theme(axis.title.x = element_text(hjust = 0))
  # labs(title = "Preset criteria (Weisman et al., 2017)")

Which capacities are attributed to which targets?

scoresplot_fun(efa_all_3, target = "all", highlight = "supernatural",
               factor_names = c("COGNITIVE", "BODILY", "SOCIAL-EMOTIONAL")) +
  scale_fill_manual(values = c("#e41a1c", "#4daf4a", "#377eb8")) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 1,
                                   color = c(rep("black", 8),
                                             rep("#984ea3", 2)),
                                   face = c(rep("plain", 8),
                                            rep("bold", 2))))

scoresplot_fun(efa_all_3, target = c("ghosts", "God", "children")) + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

itemsplot_fun(efa_all_3, target = c("ghosts", "God", "children")) + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

Comparing conceptual structures

We’ll extract 4 factors from all samples, to keep things simple.

efa_us_ad_4 <- fa(d_us_ad_pilot, nfactors = 4, rotate = "oblimin")
efa_gh_ad_4 <- fa(d_gh_ad, nfactors = 4, rotate = "oblimin")
efa_gh_ch_4 <- fa(d_gh_ch, nfactors = 4, rotate = "oblimin")
efa_th_ad_4 <- fa(d_th_ad, nfactors = 4, rotate = "oblimin")
efa_th_ch_4 <- fa(d_th_ch, nfactors = 4, rotate = "oblimin")
efa_vt_ad_4 <- fa(d_vt_ad, nfactors = 4, rotate = "oblimin")
efa_vt_ch_4 <- fa(d_vt_ch, nfactors = 4, rotate = "oblimin")

plot_us_ad_4 <- heatmap_fun(efa_us_ad_4) + 
  guides(fill = "none") + 
  labs(x = "US: adults") +
  theme(axis.title.x = element_text(hjust = 0))

plot_gh_ad_4 <- heatmap_fun(efa_gh_ad_4) + 
  guides(fill = "none") + 
  labs(x = "Ghana: adults") +
  theme(axis.title.x = element_text(hjust = 0))

plot_gh_ch_4 <- heatmap_fun(efa_gh_ch_4) + 
  guides(fill = "none") + 
  labs(x = "Ghana: children") +
  theme(axis.title.x = element_text(hjust = 0))

plot_th_ad_4 <- heatmap_fun(efa_th_ad_4) + 
  guides(fill = "none") + 
  labs(x = "Thailand: adults") +
  theme(axis.title.x = element_text(hjust = 0))

plot_th_ch_4 <- heatmap_fun(efa_th_ch_4) + 
  guides(fill = "none") + 
  labs(x = "Thailand: children") +
  theme(axis.title.x = element_text(hjust = 0))

plot_vt_ad_4 <- heatmap_fun(efa_vt_ad_4) + 
  guides(fill = "none") + 
  labs(x = "Vanuatu: adults") +
  theme(axis.title.x = element_text(hjust = 0))

plot_vt_ch_4 <- heatmap_fun(efa_vt_ch_4) + 
  guides(fill = "none") + 
  labs(x = "Vanuatu: children") +
  theme(axis.title.x = element_text(hjust = 0))

plot_grid(#plot_us_ad_4, 
          plot_gh_ad_4, plot_gh_ch_4, plot_th_ad_4, 
          plot_th_ch_4, plot_vt_ad_4, plot_vt_ch_4, 
          ncol = 2, labels = c("A", "B", "C", "D", "E", "F", "G"))

Counts

paste("US adults:", nrow(d_us_ad_pilot))
paste("GH adults:", nrow(d_gh_ad))
paste("GH children:", nrow(d_gh_ch))
paste("TH adults:", nrow(d_th_ad))
paste("TH children:", nrow(d_th_ch))
paste("VT adults:", nrow(d_vt_ad))
paste("VT children:", nrow(d_vt_ch))

paste("Non-US:", sum(nrow(d_gh_ad), nrow(d_gh_ch),
                     nrow(d_th_ad), nrow(d_th_ch),
                     nrow(d_vt_ad), nrow(d_vt_ch)))